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Abstract 

We propose a new refinement algorithm to generate size-optimal quality-guaranteed 
\ Delaunay triangulations in the plane. The algorithm takes 0(n log n + m) time, where 

rj • n is the input size and m is the output size. This is the first time-optimal Delaunay 

refinement algorithm. 

O 



1 Introduction 

> . 

Geometric domain discretizations (i.e., meshing) are essential for computer-based simula- 
tions and modeling. It is important to avoid small (and also very large) angles in such 
discretizations in order to reduce numerical and interpolation errors |SF73j . Delaunay tri- 
angulation maximizes the smallest angle among all possible triangulations of a given input 
and hence is a powerful discretization tool. Depending on the input configuration, however, 
Delaunay triangulation can have arbitrarily small angles. Thus, Delaunay refinement al- 
gorithms which iteratively insert additional points were developed to remedy this problem. 
There are other domain discretization algorithms including the quadtree-based algorithms 
|BEG94[ IMVOOj and the advancing front algorithms |Loh96j . Nevertheless, Delaunay re- 
finement method is arguably the most popular due to its theoretical guarantee and perfor- 
mance in practice. Many versions of the Delau nay refin ement is suggested in the literature 



|Che89bl lEHnTl IMITM [MPW03] |Rup93[ IShe971 |Ung04 



The first step of a Delaunay refinement algorithm is the construction of a constrained or 
conforming Delaunay triangulation of the input domain. This initial Delaunay triangulation 
is likely to have bad elements. Delaunay refinement then iteratively adds new points to the 
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(a) (b) 
Figure 1: Circumcenter vs. off-center insertion on an airfoil model. Smallest angle in both 
meshes is 32°. Delaunay refinement with circumcenters inserts 731 Steiner points and results 
in a mesh with 1430 triangles (a). On the other hand, Delaunay refinement with off-centers 
inserts 441 points and generates a mesh with 854 triangles (b). 



domain to improve the quality of the mesh and to ensure that the mesh conforms to the 
boundary of the input domain. The points inserted by the Delaunay refinement are Steiner 
points. A sequential Delaunay refinement algorithm typically adds one new vertex at each 
iteration. Each new vertex is chosen from a set of candidates — the circumcenters of bad 
triangles (to improve mesh quality) and the mid-points of input segments (to conform to 
the domain boundary). Chew [Chc89bJ showed that Delaunay refinement can be used to 
produce quality-guaranteed triangulations in two dimensions. Ruppert |Rup93| extended the 
technique for computing not only quality-guaranteed but also size-optimal triangulations. 
Later, efficient implementations |She97j . extensions to three dimensions [DBS92, Shc97|, 
generalization of input type |She97[ IMPW03j , and parallelization of the algorithm STTJ02) 
were also studied. 

Recentl y, the s econd author proposed a new insertion strategy for Delaunay refinement 
algorithm Ung04 . He introduced the so-called off-centers as an alternative to circumcen- 
ters. Off-center of a bad triangle, like circumcenter, is on the bisector of the shortest edge. 
However, for relatively skinny triangles it is closer to the shortest edge than the circumcenter 
is. It is chosen such that the triangle formed by the endpoints of the shortest edge and the 
off-center is barely of good quality. Namely, the off-center i nsertion is a more "local" opera- 
tion in the mesh than circumcenter insertion. It is shown in Ung04 that this new Delaunay 
refinement algorithm has the same theoretical guarantees as the Ruppert's refinement, and 
hence, generates quality-guaranteed size-optimal meshes. Moreover, experimental study indi- 
cates that Delaunay refinement algorithm with off-centers inserts considerably fewer Steiner 
points than the circumcenter insertion algorithms and results in smaller meshes. For in- 
stance, when the smallest angle is required to be 32°, the new algorithm inserts about 40% 
less points and outputs a mesh with about 40% less triangles (see Figure HJ). This implies 
substantial reduction not only in mesh generation time, but also in the running time of the 
application algorithm. This new off-center based Delaunay refinement algorithm is included 
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in the fifth release of the popular TrianglqJ software. Shewchuk observed (personal com- 
munication) in this new implementation, that unlike circumcenters, computing off-centers is 
numerically stable. 

Original Delaunay refinement algorithm has quadratic time complexity |Rup93|. This 
compares poorly to the time-optimal quadtree refinement algorithm of Bern et al. [BEG94 
which runs in 0(n log n + m) time, where m is t he minim um size of a good quality mesh. The 
first improvement was given by Spielman et al. STU02 as a consequence of their paralleliza- 
tion of the Delaunay refinement algorithm. Their algorithm runs in 0(m log m log 2 (L/h)) 
time (on a single processor), where L is the diameter of the domain and h is the smallest 
feature in the input. Recently, Miller Mil04j further improved this describing a new sequen- 
tial Delaunay refinement algorithm with running time 0((n log(L/h) +m)logm). In this 
paper, we present the first time optimal Delaunay refinem ent algo rithm. As Steiner points, 



we employ off-centers and generate the same output as in Ung04 . Our improvement relies 



on avoiding the potentially expensive maintenance of the entire Delaunay triangulation. In 
particular, we avoid computing very skinny Delaunay triangles, and instead we use a scaf- 
fold quadtree structure to efficiently compute, locate and insert the off-center points. Since 
the new algorithm generates the same output as the off-center based Delaunay refinement 



algorithm given by Ungor |Ung04 , it is still a Delaunay refinement algorithm. In fact, our 
algorithm implicitly computes the relevant portions of the Delaunay triangulation. 

The rest of the paper is organized as follows: In Section El we survey the necessary 
background. In Section |3J we formally define the notion of loose pairs to identify the points 
that contribute to bad triangles in a Delaunay triangulation. Next, we describe a simple 
(but not efficient) refinement algorithm based on iterative removal of loose pairs of points. 
In Section 01 we describe the new time-optimal algorithm and prove its correctness. We 
conclude with directions for future research in Sectional 



2 Background 

In two dimensions, the input domain Q is usually represented as a planar straight line graph 
(PSLG) — a proper planar drawing in which each edge is mapped to a straight line segment 
between its two endpoints |Rup93|. The segments express the boundaries of Q and the 
endpoints are the vertices of Q. The vertices and boundary segments of Q will be referred 
to as the input features. A vertex is incident to a segment if it is one of the endpoints of 
the segment. Two segments are incident if they share a common vertex. In general, if the 
domain is given as a collection of vertices only, then the boundary of its convex hull is taken 
to be the boundary of the input. 

The diametral circle of a segment is the circle whose diameter is the segment. A point is 
said to encroach a segment if it is inside the segment's diametral circle. 

Given a domain Q embedded in 1R 2 , the local feature size of each point x G 1R 2 , denoted 
by lfsn(x), is the radius of the smallest disk centered at x that touches two non-incident 
input features. This function is proven Rup93| to have the so-called Lipschitz property, i.e., 



lfsn(x) < \fsn(y) +\\xy\\, for any two points x, y E 1R 

1 http : // www- 2 . cs . emu . edu/~quake/ triangle . html 
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Figure 2: Flower of a pair of points p and q. 



In this extended abstract, we concentrate on the case where f2 is a set of points in the 
plane contained in the square [1/3, 2/3] 2 . We denote by P the current point set maintained 
by the refinement algorithm, and by £F the final point set generated. 

Let P be a point set in TR d . A simplex r formed by a subset of P points is a Delaunay 
simplex if there exists a circumsphere of r whose interior does not contain any points in P. 
This empty sphere property is often referred to as the Delaunay property. The Delaunay 
triangulation of P, denoted Del(P), is a collection of all Delaunay simplices. If the points are 
in general position, that is, if no d+2 points in P are co-spherical, then Del(P) is a simplicial 
complex. The Delaunay triangulation of a point set can be constructed in 0(n\ogn) time 
in two dimensions Ed eOlj . 

In the design and analysis of the Delaunay refinement algorithms, a common assumption 
made for the input PSLG is that the input segments do not meet at junctions with small 
angles. Ruppert |Rup93| assumed, for instance, that the smallest angle between any two 
incident input segment is at least 90°. A typical Delaunay refinement algorithm may start 
with the constrained Delaunay triangulation |Che8 9aj of the input vertices and segments or 
the Delaunay triangulation of the input vertices. In the latter case, the algorithm first splits 
the segments that are encroached by the other input features. Alternatively, for simplicity, 
we can assume that no input segment is encroached by other input features. A prepro cessing 
algorithm, which is also parallizable, to achieve this assumption is given in 

For technical reasons, as in jRup93 , we put the input Q inside a square This is 
to avoid growth of the mesh region and insertion of infinitely many Steiner points. Let 
MB = [1/3, 2/3] 2 be the minimum enclosing square of Q. The side length of B = [0, l] 2 is 
three times that of MB. We insert points on the edges of B to split each into three. This 
guarantees that no circumcenter falls outside B. We maintain this property throughout the 
algorithm execution by refining the boundary edges as necessary. 

Radius- edge ratio of a triangle is the ratio of its circumradius to the length of its shortest 
side. A triangle is considered bad if its radius-edge ratio is larger than a pre-specified constant 
P > y/2. This quality measure is equivalent to other well-known quality measures, such as 
smallest angle and aspect ratio in two dimensions |Rup93| . Consider a bad triangle, and 
observe that it must have an angle smaller or equal to a, where a = arcsin(l/2/?) 

Table Q (in the appendix) contain a summary of the notation used in this paper. 

2 In fact the reader might find it easier to read the paper, by first ignoring the boundary, (e.g., considering 
the input is a periodic point set). 
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(a) (b) 
Figure 3: The crescent of a loose pair (p, q) is shown as the shaded region (a). If the crescent 
of pq is empty of all the other vertices then the furthest point from pq inside the leaf (shown 
as c) is the off-center of pq. Otherwise, the moonstruck of a loose pair (p, q) with non-empty 
crescent is shown as r (b). Off-center in this case is the circumcenter of pqr, shown as d . 

3 Loose Pairs vs. Bad Triangles 

In the following, we use (3 to denote the user specified constant for radius-edge ratio threshold. 
Accordingly, a denotes the threshold for small angles. 

Definition 3.1 For a pair of vertices p and q, let di(p,q) be the disk with center u such 
that pqu is a left_turn and ||wp|| = ||wg|| = Similarly, let d r (p, q) be the disk with 

center v such that pqv is a right_turn and \\vp\\ = \\vq\\ = ft\\pq\\- We call the union of the 
disks d[(p,q) and d r (p, q) the flower of pq. Moreover, d/(p, q) is called the left leaf of the 
flower and d r (p, q) is called the right leaf of the flower. 

A pair of vertices (p, q) in P is a loose pair if either the left or the right leaf of the flower 
of pq is empty of vertices. 

Let (p, q) be a loose pair due to an empty left (resp. right) leaf, and c be the furthest 
point from pq on the boundary of the leaf. See Figure El (a). Let d be the disk centered at 
c having p and q on its boundary. We call the region d \ di(p, q) (resp. d \ d r (p, q)) the left 
(resp. right) crescent of pq and denote it by crescent; (pq) (resp. crescent r (pg)). 

Crescent of a loose pair (p, q) may or may not be empty of all the other vertices. In the 
latter case, the moonstruck of pq is the vertex r inside the crescent such that the circumdisk 
of pqr is empty of all the other vertices, see Figure El (b). 

Lemma 3.2 There exists a loose pair in a point set P if and only if the minimum angle in 
the Delaunay triangulation of P is smaller than or equal to a. 

Proof: If there exists a loose pair (p, q) then pq is a Delaunay edge. Moreover the triangle 
incident to edge pq on the side of the empty leaf must be bad, with an angle smaller than a. 

For the other direction, let pqr be the bad triangle with the shortest edge. Without loss 
of generality assume pqr is a right_turn. If the right leaf of pq is empty then (p, q) is a loose 
pair and we are done. Otherwise, let s be the first point a morphing from the circumsphere 
of pqr to right leaf of pq hits (fixing the points p and q). Then both ps and sq are shorter 
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Algorithm 1 Loose Pair Removal 
Require: A point set f2 in R 2 and (5 

Ensure: A Steiner triangulation of Q where all triangles have radius-edge ratio at most /3 
Let P = Q 

while there exists a loose pair (p, q) in P do 

Insert the off-center of pq into P 
end while 

Compute and Output the Delaunay triangulation of the resulting point set 



features then pq and are incident to a bad Delaunay triangle. This is a contradiction to the 
minimality of the \ \pq\\- ■ 
This lemma suggest the refinement method depicted in Algorithm ^ Note that each 
loose pair corresponds to a bad triangle in the Delaunay triangulation of the growing point 
set. Hence, this algorithm is sim ply ano ther way of stating the Delaunay refinement with 



off-centers algorithm presented in Ung04 . Ungor showed that the Delaunay refinement with 
off-centers algorithm terminates and the resulting point set is size-optimal. In order to give 
optimal time bounds we will refine this algo rithm in the next section. 

They state that during the refinement 



The next two lemmas follow directly from Ung04 



process we never introduce new features that are smaller than the current loose pair being 
handled. 



Lemma 3.3 Let P be a point-set, (p, q) be a loose pair of P, and P' be the set resulting 
for inserting the off-center of pq. We have for any x e P, that if\fsp>(x) < lfsp(x), then 
lfsp'(x) > ||pg|| . 

Lemma 3.4 Let Q be the input point set, and let P be the current point set maintained by 
the refinement algorithm depicted in Algorithm{][ Let 5" denote the point set generated by 
Algorithm^ Then for any point p in the plane, we have throughout the algorithm execution 
that lfen(p) > lfsp(p) > lfss-(p) > c shrink \fsn(p) , where c shrink > is a constant. 

Lemma 13.31 suggests a natural algorithm for generating a good cloud of points. Since 
inserting a new off-center can not decrease the smallest feature of the point cloud, it is natural 
to first handle the shortest loose pair first. Namely, repeatedly find the smallest loose pair, 
insert its off-center, till there are no loose pairs left. Because the domain is compact, by a 
simple packing argument it follows that this algorithm terminates and generates an optimal 
mesh. This is one possible implementation of (the generic) Algorithm ^ 

Implementing this in the naive way, is not going to be efficient. Indeed, first we need 
to maintain a heap sorted by the lengths of the loose pairs, which is already too expensive. 
More importantly, checking if a pair is loose requires performing local queries on the geometry 
which might be too expensive to perform. 

We will overcome these two challenges by handling the loose pairs using a weak ordering 
on the pairs. This would be facilitated by using a quadtree for answering the range search- 
ing queries needed for the loose pairs determination. In particular, our new algorithm is 
just going to be one possible implementation of Algorithm Q and as such Lemma 13.31 and 
Lemma f3.4l hold for it. 
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Algorithm 2 : Delaunay .Refinement 
Require: A point set Q C [1/3, 2/3] 2 G R 2 and (3 

Ensure: A Steiner triangulation of Q where all triangles have radius-edge ratio at most /3 
Split each edge of the square [0, l] 2 into three segments. 
Construct a balanced quadtree QT of fl using [0, l] 2 for the root node. 
Insert all the nodes of QT into a heap TCP, sorted from smallest node to largest. 
Initialize all vertices to be active. 
Let prevJ be the depth of the smallest QT node, 
while TiV is not empty do 

□ «- extractMin(W). 

i = depth(D). 

if % < prevJ then 

Move all vertices in level prevJ which are also active in level i to the ith level. 
prevJ = i. 
end if 

for every active vertex p G P fl □ do 

for every active vertex q G P such that \\pq\\ < c reac hsize(n) do 
if pq is loose then 

Insert the r = off-center (p, q), and store the r in QT in a cell □' as low as 
possible, such that q ou ,size(D') <||pr|| < c up size(D') and size(D') > size(D) 
end if 

for every node in the same level of □ that had a point inserted into it, because of 
the above step, reinsert it into the heap TCV. 
end for 
end for 
end while 

Compute and Output the Delaunay triangulation of the resulting point set 



4 Efficient Algorithm Using a Quadtree 

We construct a balanced quadtree for Q, using the unit square as the root of the quadtree. In 
the following, P denotes the current point set, as it grows during the algorithm execution. Let 
be the final point set generated. A balanced quadtree has the property that two adjacent 
leaves have the same size up to factor two. A balanced quadtree can be constructed, in 
0(n log n+m) time, where m is the size of the output, see BEG94]. Such a balanced quadtree 
also approximates the local feature size of the input, and its output size m is (asymptotically) 
the size of the cloud of points we need to generate. In the constructed quadtree we maintain, 
for each node, pointers to its neighbors in its own level, in the quadtree, and to its neighbors 
in the levels immediately adjacent to it. 

The new algorithm is depicted in Algorithm |21 For the time being, consider all points 
to be active throughout the execution of the algorithm. Later, we will demonstrate that it 
is enough to maintain only very few active points inside each cell, thus resulting in a fast 
implementation. We show that each quadtree node is rescheduled into the heap at most a 
constant number of times, implying that the algorithm terminates. 



7 



In Algorithm EJ collecting the active points around a cell □, checking whether a pair is 
active, or finding the moonstruck point of a pair is done by traversing the cells adjacent to 
the current cell, using the boundary pointers of the well-balanced quadtree. We will show 
that all those operations takes constant time per cell. 

One technicality that is omitted from the description of Algorithm |21 is that we refine 
the boundary edges of the unit square by splitting such an edge in the middle, if it is being 
encroached upon. Since the local feature size on the boundary of the unit square is 0(1), by 
Lemma (3.41 As such, this can automatically handled every time we introduce a new point, 
and it would require 0(1) time for each insertion. This guarantees that no point would be 
inserted outside the unit square. Note, that in such a case the encroaching new vertex is not 
being inserted into the point set (although it might be inserted at some later iteration). 

4.1 Proof of Correctness 

The proof of correctness is by induction over the depth of the nodes being handled. We 
use dgj to denote the depth of the quadtree QT. In the kth stage of the execution of the 
algorithm, it handles all nodes of depth (d,Qq — k) in the tree. Next, the algorithm handles 
all nodes of depth (g?qt — (k + 1)), and so on. 

By the balanced quadtree construction [BEG94J, for every leaf □ of the quadtree QT, 
and every point p G P R □, we have Q OW) size(D) < lfsp(p) < c up size(D), where ci ow and c up 
are prespecified constants such that c up > 1c\ ow . In particular, the value of c up and c\ ow is 
determined by the initially constructed quadtree. 

Lemma 4.1 Let P be the current point set maintained by Algorithm^ and let r be an off- 
center of a loose pair (p, q) in P. Let be the quadtree node that the point r is inserted into. 
We have C[ ow ■ size(D') < lfsp(r) < c up ■ size(D'). In particular, c' low ■ size(D') < lfs^(r) < 
Cup ' size(d ) ; where c^ ow = c s hrink ■ ciow 

Proof: The claim follows from the explicit condition used in the insertion part of the 
algorithm. Observe that since lfsp(r) = ||pr|| > lfsp(p) > C[ ow ■ size(D), where □ is the cell of 
the quadtree containing the point p. As such, a node □' that contains r and is in the same 
level as □, will have Q ow -size(D') < lfsp(p) < lfsp(r), by induction. If lfsp(r) < c up -size(D / ) 
then we are done. Otherwise, lfsp(r) > c up ■ size(D') implying that ci ow ■ size (parent (□')) < 
Cup ' (□') < lfsp(r), since c up > 2q ow . Thus, set □' <— parent (□') and observe that 
Q ow -size(D') < lfsp(r), as such we can continue climbing up the quadtree till both inequalities 
hold simultaneously. 

The second part follows immediately from Lemma f3 .41 ■ 
Lemma 4.2 Off-center insertion takes 0(1) time. 

Proof: Let r be an off-center of a loose pair (p,q), and let □ and □' be the cells of the 
quadtree containing p and r, respectively. By Lemma 14.11 lfsp(p) = 0(size(D)). By the 
algorithm definition, we have = ©(size(D)). As such, ||pr|| = 0(||pg||) = ©(size(D)). 
Namely, lfsp(r) = ||pr|| = 0(size(d)). Namely, in the grid resolution of □, the points p 
and r are constant number of cells away form each other (although r might be stored a 
constant number of levels above □). Since every node in the quadtree have cross pointers 
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to its immediate neighbors in its level (or the above level), by the well-balanced property of 
the quadtree. It follows that we can traverse from □ to □' using constant time. ■ 

Lemma 4.3 The shortest loose pair of P in the beginning of the ith stage is of length at 
least r ]i = c' low /2 d ^- i + 1 . 

Proof: For % = 1, the claim trivially holds by the construction of the balanced quadtree 
using [BEC^94] . Now, assume that the claim holds for i — 1, . . . , k. We next show that the 
claim holds for i — k + 1. Specifically, it holds at the end of the kth stage. 

Suppose for the sake of contradiction that there exists a loose pair (p, q) shorter than r/k- 
Assume, without lose of generality, that p was created after q by the algorithm, and let □ 
be the node of QT that contains p. If the depth of □ is dgq — k, then since \\pq\\ < r]k < 
Creachsize(n), we have that the algorithm handled the point p in □, it also considered the 
pair (p, q) and inserted its off-center. So, the pair (p, q) is not loose at the end of the fcth 
stage. 

If the depth of □ is larger then dgq — k then lfs(p) was larger than Q ou; size(n) when p 
was inserted. As such, lfsp(p) > c' low size(D) > rjk, which is a contradiction, since any loose 
pair that p participates in must be of length at least lfsp(p). ■ 

Note that when the algorithm handles the root node in the last iteration, it "deteriorates" 
into being Algorithm ^ executed on the whole point set. Hence Lemma 14 .HI implies the 
following claim. 

Claim 4.4 In the end of the execution of Algorithm there are no loose pairs left in jF, 
where 3 is the point set generated by the algorithm. 

4.2 How the refinement evolves 

Our next task, is to understand how the refinement takes place around a point, and form a 
"protection" area around it. In particular, the region around a point p G P with lfsn(p) is 
going to be effected (i.e., points inserted into it), starting when the algorithm handles cells 
of level i, where 1/2* « lfso(p). Namely, the region around p might be refined in the next 
few levels. However, after a constant number of such levels, the point p is surrounded by 
other points, and p is not loose with any of those points. As such, the point p is no longer 
a candidate to be in a loose pair. To capture this intuition, we prove that this encirclement 
process indeed takes place. 

We define the gap of a vertex x G P (which is not a boundary vertex), denoted by gap(x), 
as the ratio between the radius of the largest disk that touches x and does not contain any 
vertex inside, and the lfsp(x). 

Lemma 4.5 For a vertex w G P, if gap (w) > c 9 = (2/3) 7r / Q+1 , then there exists a loose pair 
of P of length < c gbu ■ \fs(w), where c gbu = (2/?) 7r/af . 

Proof: Assume that w is strictly inside the bounding square B. Proof for the case where 
w is on the boundary of B is similar and hence omitted. Let T 1; T 2 , . . . , T m be the Delaunay 
triangles incident to w and Ui,u 2 , ■ ■ ■ ,u m be the Delaunay neighbors of w. Note that if 
||wuj|| > 2 / 9||wMj_ 1 ||, then the left flower of wu^i must be empty and hence wu^i is a 
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loose pair. Similarly, if ZwUiU^i < a then (w,iti_i) is a loose pair. If Zwu^Ui-x > a and 
Z.UiWUi^i < a then by the law of sines, it must be that < ||«j_iit;||, and is 

facing an angle smaller than a, and as such it is a loose pair. 

Suppose that none of the triangles T%, . . . , T m have a loose pair on their boundary. Then, 
it must be that m < 2n/a, since the angle Z.UiWUi + \ > a, for i = l,...,m. But then, 
IKHI < (2/3) i ~ 1 \fs(w) andllwi^H < (2/3) ro -* +1 lfs(>). It follows that \\ Ui w\\ < (2/3) w / a lfe(u>), 
for i — 1, . . . , m. Since, all the angles in Tj are larger than a, it follows that circumcircle of 
Tj is of radius < < /3(2/3) 7r//a lfs(w). But then, the gap around w, is at most P(2f3) 7T / a . 

A contradiction, since gap(iu) = c g > (3{2(3Y^ a . 

Thus, one 7\, . . . , T m must be bad. Arguing as above, one can show that the first such 
triangle, has a loose pair of length < (2/?) 7r//Q lfs(w), as claimed. ■ 

Lemma [4.51 implies that if we handle all loose pairs of length smaller than £, then all the 
points having a big gap, must be with local feature size Q(£). 

Lemma 4.6 Let P be a point set such that all the loose pairs are of length at least tjc\, for 
a constant c\>2. Let (p, q) be a loose pair of length i with a non-empty crescent, and let w 
be its moonstruck neighbor. Then, lfs(w) > eic b . 

Proof: Since w is a moonstruck point of a loose pair of length at least i there is an 
empty ball of radius at least f3l touching w. Hence, gap(w) > ^^y- We consider two 

cases. If jjr^y > c g , then by Lemma l4~H| there exist a loose pair of size at most c g b U lfs(w). 
However, all loose pairs are of length > £/c\, and it follows that c g b u ^k(w) > ijc\. Hence, 
lfs(w) > £/(ciC g b u ). On the other hand, if < c g then, 

lfs(w) > — > — > > , 

Cg Cg Cgbu C\Cgfy U 

since c g b u > c g and (3 >1. ■ 



4.3 Managing Active Points 

As we progress with execution of the algorithm, the results of the previous section imply 
that a vertex with relatively small feature size cannot participate in a loose pair, nor be a 
moonstruck point. So, in the evolving quadtree we do not maintain such set of points that 
play no role in the later stages of the algorithm execution. This facilitates an efficient search 
for finding the loose pairs and moonstruck points as shown in the rest of this section. For 
each vertex, size of its insertion cell gives a good approximation of its feature size. We use 
this to determine the lifetime of each vertex in our evolving quadtree data structure. 

The activation depth of an input vertex p, denoted by p, is the level of the initial quadtree 
leaf containing p. For a Steiner point p, the activation depth is the depth of the cell p is 
inserted into. 

Lemma 4.7 A vertex p can not be a loose pair end or a moonstruck point, when the algo- 
rithm handles level of depth <p — c span; where c span = lg Cgb / " c " p + 1. 
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Proof: When the point p was created, we had lfsp(p) < c up /2^, by Lemma 14.11 Now, 
if p is an endpoint of a loose pair in depth m < p in the quadtree, it must be that the 
length t of this pair is at least c' low /2 m , by Lemma IP1 Since the local feature size lfsp is a 
non-increasing function as our algorithm progresses, it follows that 

o (J /o m _ (J 
lfsp(p) c up /2P c up 

If gap(p) < c g then 2^~ m ■ ^» < c g . Implying that p - m < lg ^ s -. 

By Lemma T4.5[ if gap(p) at any point in the algorithm becomes larger than c g , then there 
exists a loose pair of length < c g b u • lfsp(p). But all such pairs are handled in level > t in the 
quadtree, where c g b u • lfsp(p) > c' low /2 t+1 by Lemma ED Thus, c gbu c up /2 p > c gbu Us P (p) > 
c' low /2 t+1 . Implying that 

t > p=p-\g y f - 1. 

C low 

This implies, that when the algorithm handles cells of depth 1, . . . , p, we have that the vertex 
p can not participate directly in a loose pair. 

If p is not loose pair end, but is a moonstruck point for a loose pair, then 

° U P \ lf„ l„\ \ C low 



> lfep(p) > 



2P _ _ 2t+lc ^ u 

by Lemma [4.61 and Lemma [4.31 This in turn implies that t > p — lg Cg6 " 6 " p — 1. ■ 

low 

Definition 4.8 A point p is active at depth i, if p > z > p — c span , where c span is a constant 
specified in Lemma f4. 71 

Note, that the algorithm can easily maintain the set of the active points. Lemma 14.71 
implies that only active points are needed to be considered in the loose pair computation. 

Observation 4.9 During the off-center insertion any new loose pairs introduced are at least 
the size of the existing loose pairs. 

Lemma 4.10 At any stage i, the number of active points inside a cell at level dQq — i is a 
constant. 

Proof: This is trivially true in the beginning of the execution of the algorithm, as the 
initial balanced quadtree has at most a constant number of vertices in each leaf. Later on, 
Lemma f4. 71 implies that when a point p is being created, with I = lfs(p), then its final local 
feature size is going to be Q(£). To see that, observe that when p was created, the algorithm 
handled loose pairs of size Q(i). From this point on, the algorithm only handle loose pairs 
that are longer (or slightly shorter) than I. Such a loose pair, can not decrease the local 
feature size to be much smaller than £, by Lemma 13.31 

This implies that when p is being created, we can place around it a ball of radius f2(lfs(j»)) 
which would contain only p in the final generated point set. Since p becomes inactive c span 
levels above the level it is being created, it follows that a call in the quadtree can contain at 
most a constant number of such protecting balls, by a simple packing argument. ■ 
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4.4 Efficient Implementation Details 

The above discussion implies that during the algorithm execution, we can maintain for every 
quadtree node a list of constant size that contains all the active vertices inside it. When 
processing a node, we need to extract all the active points close to this cell □. This requires 
collecting all the cells in this level, which are constant number of cells away from □ in this 
grid resolution. In fact, the algorithm would do this point collection also in a constant 
number of levels above the current level, so that it collects all the Steiner points that might 
have been inserted. Since throughout the execution of the algorithm we maintain a balanced 
quadtree, we have from every node, pointers to its neighbors either in its level, or at most 
one level up. As such, we can collect all the neighbors of □ in constant distance from it 
in the quadtree, in constant time, and furthermore, extract their active points in constant 
time. Hence, handling a node in the main loop of Algorithm El takes constant time. 

We need also to implement the heap used by the algorithm. We store nodes in the 
heap TiV, and it extracts them according to their depth in the quadtree. As such, we can 
implement it by having a separate heap for each level of the quadtree. Note that the local 
feature size of a vertex when inserted into a quadtree node is within a constant factor of the 
size of the node. Hence a node can be rescheduled in the heap at most a constant number 
of times. For each level, the heap is implemented by using a linked list and a hash-table. 
Thus, every heap operation takes constant time. 

4.5 Connecting the Dots 

We shall also address how to perform the final step of Algorithm that is computing the 
Delaunay triangulation of the resulting point set £F. This can be done by re-executing a 
variant of the main loop of Algorithm|2]on 5", which instead of refining the point set, reports 
the Delaunay triangles. We use a similar deactivation scheme to ignore vertices whose all 
Delaunay triangles are reported. Since 3 is a well-spaced point set, for a pair of nearby 
active vertices we can efficiently compute whether the two makes a Delaunay edge and if so 
locate also the third point that would make the Delaunay triangle. It is straightforward but 
tedious to argue that the running time of this algorithm is going to be proportional to the 
running time of Algorithm El 

4.6 Analysis 

The initial balanced quadtree construction takes 0(n\ogn + m) time, where m is the size 
of the resulting quadtree [BEG94J. This quadtree has the property that the size length of 
a leaf is proportional to the local feature size. This in turn implies the value of ci ow and 
c up , which in turn guarantees that no new leafs would be added to the quadtree during the 
refinement process. 

Furthermore, the point set generated by Algorithm El has the property that its density is 
proportional to the local feature size of the input. Namely, the size of the generated point-set 
is 0(m). Since all the operations inside the loop of Algorithm Stakes constant time, we can 
charge them to either the newly created points, or to the relevant nodes in the quadtree. 
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This immediately implies that once the quadtree is constructed, the running time of the 
algorithm is 0(m). 

Theorem 4.11 Given a set Q of n points in the plane the Delaunay refinement algorithm 
(depicted in Algorithm^ computes a quality- guaranteed size-optimal Steiner triangulation 
ofQ, in optimal time 0(nlogn + m) , where m is the size of the resulting triangulation. 



5 Conclusions 

We presented a time-optimal algorithm for Delaunay refinement in the plane. It is important 
to note that the output of this new algorithm is the same as that of the off-center based 



Delaunay refinement algorithm given in Ung04 , which outperforms the circumcenter based 



refinement algorithms in practice. The natural open question for further research is extending 
the algorithm in three (and higher) dimensions. We believe that extending our algorithm to 
handle PSLG in the plane is doable (with the same time bounds), but is not trivial, and it 
would be included in the full version of this paper. 

We note that when building the initial quadtree, we do not have to perform as many re- 
finement steps as used in the standard quadtree refinement algorithm of Bern et al. [BEG94J. 
While their algorithm considers a quadtree cell with two input vertices crowded and splits 
it into four, we are perfectly satisfied with a balanced quadtree as long as the quadtree ap- 
proximates the local feature size within a constant and hence the number of features in a cell 
is bounded by a constant. This difference in the depth of the quadtree should be exploited 
for an efficient implementation of our algorithm (this effects the values of the constants ci ow 
and c up ). 

Parallelization of quadtree based methods are well understood BET99J, while design of 
a theoretically optimal and pra ctical parallel Delaunay refinement algorithm is an ongoing 



research topic |STU02[ lSTU04j . We believe our approach of combining the strengths of 
quadtrees as a domain decomposition scheme and Delaunay refinement with off-centers will 
lead to a good parallel solution for the meshing problem. 
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Notation 


Value 


Comment 


n 




Input point set. 


p 




Current point set maintained by the algorithm. 






Final point set 


P 


> V2 


radius-edge ratio threshold for bad triangles 


a 


arcsin( 1/2/3) 


small angle threshold bad triangles 


C 9 
Cgbu 
Clow 
Cup 


(2f3Y' a+1 
(2/3)f/« 

— 2ci ow 


Vertex with larger gap than c g participates in a loose pair. 
Blowup of lfs for a lose pair around a vertex with large gap. 
Lower bound on the lfs of a point inside a leaf of the initial quadtree. 
Upper bound on the lfs of a point inside a leaf of the quadtree. 


Crcach 


2,C U pCgl) U 


Crcachsize(D) is an upper bound on the length of a 
(relevant) lose pair involving a point p stored in a cell 
□. For correctness, it required that c reac h > 2\/2. 


Cshrink 


> 


For any point x, we have lfs^(x) > c s h r ink • lfen(^)- 


r > 

How 


Clow ' Cshrink 


Lower bound on the lfs of a point p stored inside a 
node □ of the quadtree through the algorithm execu- 
tion. Namely, lfsgr(p) > c' low size(n) . 



Table 1: Notation used in the paper. 
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